#Replication code, `Banditry or Business? Rebel Labor Markets and State Economic Intervention'
#C. Estancona
#08/06/2021`

rm(list=ls(all=TRUE))
library(foreign)
library(lme4)
library(readstata13)
library(ggplot2)
library(stargazer)
library(xtable)
library(miceadds)
library(DataCombine)
library(dplyr)


#read in data
Colombia_1<-read.csv("Colombia_credits.csv")

#Create DV, log of agricultural credits
Colombia_1$ln_ba_peq_vr<-log(Colombia_1$ba_peq_vr+1)

#Create dummy for if credits are provided or not
Colombia_1$d_ba<-ifelse(Colombia_1$ba_tot_vr>0, 1,0)

#Create dummy for coca presence
Colombia_1$coca_pos<-ifelse(Colombia_1$D_Coca ==1 | Colombia_1$H_coca>0, 1, 0)

#create dummy by municipality for FARC possibility (for summaries)
Colombia_1<-Colombia_1 %>% 
  group_by(codmpio) %>% 
  mutate(FARC_pres=max(FARC, na.rm=T))

#necessary to use with slide() below
Colombia_1<-as.data.frame(Colombia_1)


#to lead credits (credits predicted by all other variables in year prior)
Colombia_slid <- slide(Colombia_1, Var = "ln_ba_tot_vr", GroupVar = "codmpio",
                   slideBy = 1)

Colombia_slid_2 <- slide(Colombia_1, Var = "ln_ba_peq_vr", GroupVar = "codmpio",
                       slideBy = 1)

Colombia_slid_4 <- slide(Colombia_1, Var = "d_ba", GroupVar = "codmpio",
                         slideBy = 1)

####################################
#######Table 1 in Manuscript 

Colombia_2<-Colombia_1[which(Colombia_1$ano>=2000),]

#average value by municipality characteristics
xtable(with(Colombia_2, tapply(ba_tot_vr, list(FARC=FARC_pres,coca_pos=coca_pos), mean)))

##########################
######Main table in text (table 2)

#dummy, coca presence (M1)
de_lag<-lm(ln_ba_tot_vr~FARC*coca_pos+as.factor(codmpio)+as.factor(ano), data=Colombia_slid)

#continuous coca measure (M2)
de_lag_2<-lm(ln_ba_tot_vr~FARC*ln_H_coca+as.factor(codmpio)+as.factor(ano), data=Colombia_slid)

#to TeX
stargazer(de_lag, de_lag_2, omit=c("codmpio", "ano"))

############################
###################Appendix Tables
############################

m8.2<-lm(ln_ba_tot_vr1~FARC+FARC*D_Coca+D_Coca+AUC+g_terreno, data=Colombia_slid)
summary(m8.2) #Appendix Table A1

m8<-lm(ln_ba_tot_vr1~FARC+FARC*ln_H_coca+ln_H_coca+AUC+g_terreno, data=Colombia_slid) 
summary(m8) #Appendix Table A2, Model 1

m8.1<-lm(ln_ba_peq_vr1~FARC+FARC*+ln_H_coca+ln_H_coca+AUC+g_terreno, data=Colombia_slid_2)
summary(m8.1) #Appendix Table A2, Model 2

m8.3<-lm(ln_ba_tot_vr1~FARC+FARC*+ln_H_coca+ln_H_coca+AUC+g_terreno+disbogota, data=Colombia_slid)
summary(m8.3) #Appendix Table A2, Model 3

##############Logit models in appendix (Table A3)

m8.4<-glm(d_ba1~FARC+FARC*+ln_H_coca+ln_H_coca+AUC+g_terreno, data=Colombia_slid_4, family="binomial")
summary(m8.4)

m8.5<-glm(d_ba1~FARC+FARC*+D_Coca+D_Coca+AUC+g_terreno, data=Colombia_slid_4, family="binomial")
summary(m8.5)

#models to TeX
stargazer(m8, m8.1, m8.3)
stargazer(m8.4, m8.5)

############Fixed effect model - Appendix A4
m8.6<-lm(ln_ba_tot_vr1~FARC+FARC*+ln_H_coca+ln_H_coca+AUC+g_terreno+as.factor(ano), data=Colombia_slid)
summary(m8.6)  
#model to TeX
stargazer(m8.6)

####################Appendix Table A5:
#create a single year-department indicator for clustered standard errors 
Colombia_slid$department_year <- as.factor(paste(Colombia_slid$coddepto, Colombia_slid$ano))
Colombia_slid$department_year

#try main model with this as a clustered se
cluster_mod_2<-lm(ln_ba_tot_vr~FARC*coca_pos+as.factor(codmpio)+as.factor(ano), data=Colombia_slid)
summary(cluster_mod_2, cluster="department_year") 

#to TeX
stargazer(cluster_mod_2, se=list(coef(summary(cluster_mod_2,cluster = c("department_year")))[, 2]), omit=c("codmpio", "ano")) 



#########################################
##############Plots and summaries

############################
#######Point Estimate plot - Figure 1, Main Text

x_vec_coca_2<-rep(c(0,1))
x_vec_FARC<-rep(c(0,1))

new_dat_m8.2<-data.frame(expand.grid(x_vec_FARC, x_vec_coca_2, 1, mean(na.omit(Colombia_slid$g_terreno))))

colnames(new_dat_m8.2)<-c("FARC", "D_Coca", "AUC", "g_terreno")


fit_m8.2<-predict(m8.2, newdata=new_dat_m8.2, interval="confidence")

plot_dat_m8.2<-data.frame(fit_m8.2, new_dat_m8.2, c(1,2,3,4))
colnames(plot_dat_m8.2)[8]<-"index"
plot_dat_m8.2_nf<-plot_dat_m8.2[which(plot_dat_m8.2$FARC==0),]
plot_dat_m8.2_f<-plot_dat_m8.2[which(plot_dat_m8.2$FARC==1),]

plot_8.2_ind<-ggplot(plot_dat_m8.2, aes(x=as.factor(index), y=exp(fit), colour=as.factor(D_Coca))) + 
  geom_point(aes(x=as.factor(index), y=exp(fit), size=0.2), show.legend=FALSE)+
  geom_errorbar(aes(ymin=exp(lwr), ymax=exp(upr)))+
  scale_color_manual(values=c("grey58", "black"), name="Coca Presence", labels=c("No", "Yes"))+
  ylim(0,400)+
  scale_x_discrete(name="Municipality Type", labels=c( "No FARC", "FARC", "No FARC", "FARC"))+
  ylab("Value of Agricultural Credits")+
  theme_bw()+
  theme(axis.text.x=element_text(angle=45, hjust=1))

#######################################

######Appendix plot 6
#########
x_vec_coca<-seq(min(na.omit(Colombia_slid$ln_H_coca)), max(na.omit(Colombia_slid$ln_H_coca)), length.out=5000)

new_dat_m8_0<-data.frame(0, x_vec_coca, 1, mean(na.omit(Colombia_slid$g_terreno)))
new_dat_m8_1<-data.frame(1, x_vec_coca, 1, mean(na.omit(Colombia_slid$g_terreno)))

colnames(new_dat_m8_0)<-c("FARC", "ln_H_coca", "AUC", "g_terreno")
colnames(new_dat_m8_1)<-c("FARC", "ln_H_coca", "AUC", "g_terreno")


fit_m8_1.0<-predict(m8, newdata=new_dat_m8_0, interval="confidence")
fit_m8_1.1<-predict(m8, newdata=new_dat_m8_1, interval="confidence")

plot_dat_m8<-data.frame(fit_m8_1.0, fit_m8_1.1, new_dat_m8_1)

ggplot(data=plot_dat_m8)+ 
  labs(x="ln(Hectares of Coca)", y="ln(Value of Agricultural Credits Provided)")+
  geom_line(aes(x=ln_H_coca, y=fit_m8_1.0[,1]),linetype="solid", color="black")+
  geom_line(aes(x=ln_H_coca, y=fit_m8_1.1[,1]), linetype="longdash", color="grey58")+
  geom_ribbon(aes(x=ln_H_coca, ymin=fit_m8_1.0[,2], ymax=fit_m8_1.0[,3], fill="Non-FARC Municipalities"), alpha=.3)+
  geom_ribbon(aes(x=ln_H_coca, ymin=fit_m8_1.1[,2], ymax=fit_m8_1.1[,3], fill="FARC Municipalities"), alpha=.3)+
  theme(legend.title=element_text(size=14), legend.position="bottom", panel.grid.minor.x=element_blank(),  panel.grid.minor.y=element_blank()) +
  scale_fill_manual(values=c("black","grey58"), name="")+
  scale_color_manual(values=c("grey58", "black"), name="")+
  theme_bw(base_size=20)+
  theme(legend.position="bottom")
dev.off()

############################


########################
#######Descriptive statistics, continuous variables (Appendix Table A6)

descriptive_1<-with(Colombia_slid, as.data.frame(rbind(summary(log(H_coca+1)), summary(ln_ba_tot_vr), summary(g_terreno))))
descriptive_1<-descriptive_1[,1:6]
descriptive_1
rownames(descriptive_1)<-c("Ln(Hectares of Coca)", "Ln(Agricultural Credit Value)", "Gini (Land)")
xtable(descriptive_1)


############
##### Histogram, AUC presence (Appendix Figure A2)

ggplot(data=Colombia_slid, aes(x=factor(AUC))) +
  geom_bar(stat="count")+
  theme_bw(base_size=20)+
  labs(x="AUC Presence", y="Count")


############
##### Histogram, FARC presence (Appendix Figure A3)


ggplot(data=Colombia_slid, aes(x=factor(FARC))) +
  geom_bar(stat="count")+
  theme_bw(base_size=20)+
  labs(x="FARC Presence", y="Count")


